10  Centroids of transport zones

In this section we will calculate the geometric and the weighted centroids of transport zones.

10.1 Geometric centroids

Taking the Municipalities_geo data from the previous section, we will calculate the geometric centroids, using the st_centroid() function.

Code 
library(dplyr)
library(sf)
library(mapview)

Municipalities_geo = st_read("../data/Municipalities_geo.gpkg", quiet = TRUE)

Centroid_geo = st_centroid(Municipalities_geo)

This creates points at the geometric center of each polygon.

Code 
mapview(Centroid_geo)
Code 
mapview(Centroid_geo) + mapview(Municipalities_geo, alpha.regions = 0) # both maps, with full transparency in polygons

But… is this the best way to represent the center of a transport zone?

These results may be biased by the shape of the polygons, and not represent where activities are. Example: lakes, forests, etc.

To overcome this, we can use weighted centroids.

10.2 Weighted centroids

We will weight centroids of transport zones by population and by number of buildings.

For this, we will need the census data (INE 2022).

Code 
Census = st_read("../data/census.gpkg", quiet = TRUE)

mapview(Census |> filter(Municipality == "Lisboa"), zcol = "Population")

It was not that easy to estimate weighted centroids with R, as it is with GIS software. But there is this new package centr that can help us (Zomorrodi 2024).

We need to specify the group we want to calculate the mean centroids, and the weight variable we want to use.

Code 
library(centr)
Centroid_pop = Census |> 
  mean_center(group = "Municipality", weight = "Population")

We can do the same for buildings.

Code 
Centroid_build = Census |> 
  mean_center(group = "Municipality", weight = "Buildings")

10.3 Compare centroids in a map

10.3.1 Interactive map

Code 
mapview(Centroid_geo, col.region = "blue") +
  mapview(Centroid_pop, col.region = "red") +
  mapview(Centroid_build, col.region = "black") +
  mapview(Municipalities_geo, alpha.regions = 0) # polygon limits

See how the building, population and geometric centroids of Sintra are apart, from closer to Tagus, to the rural area.

10.3.2 Static map

To produce the same map, using only plot() and st_geometry(), we need to make sure that the geometries have the same crs.

Code 
st_crs(Centroid_geo) # 4326 WGS84
st_crs(Centroid_pop) # 3763 Portugal TM06

So, we need to transform the geometries to the same crs.

Code 
Centroid_pop = st_transform(Centroid_pop, crs = 4326)
Centroid_build = st_transform(Centroid_build, crs = 4326)

Now, to use plot() we incrementally add layers to the plot.

Code 
# Plot the Municipalities_geo polygons first (with no fill)
plot(st_geometry(Municipalities_geo), col = NA, border = "black")

# Add the Centroids_geo points in blue
plot(st_geometry(Centroid_geo), col = "blue", pch = 16, add = TRUE) # add!

# Add the Centroid_pop points in red
plot(st_geometry(Centroid_pop), col = "red", pch = 16, add = TRUE)

# Add the Centroid_build points in black
plot(st_geometry(Centroid_build), col = "black", pch = 16, add = TRUE)

Static map of different centroids of Municipalities

In the next section we will use these centroids to calculate the desire lines between them.

---
format:
  pdf:
    prefer-html: true
---

# Centroids of transport zones

In this section we will calculate the geometric and the weighted centroids of transport zones.

## Geometric centroids

Taking the `Municipalities_geo` data from the previous section, we will calculate the geometric centroids, using the `st_centroid()` function.

```{r geocentroid}
#| message: false
#| warning: false
library(dplyr)
library(sf)
library(mapview)

Municipalities_geo = st_read("../data/Municipalities_geo.gpkg", quiet = TRUE)

Centroid_geo = st_centroid(Municipalities_geo)
```

This creates points at the geometric center of each polygon.

quarto-executable-code-5450563D

```r
#| message: false
#| fig-format: png
mapview(Centroid_geo)
mapview(Centroid_geo) + mapview(Municipalities_geo, alpha.regions = 0) # both maps, with full transparency in polygons
```

But...
is this the best way to represent the center of a transport zone?

These results may be biased by the shape of the polygons, and not represent where activities are.
Example: lakes, forests, etc.

To overcome this, we can use **weighted centroids**.

## Weighted centroids

We will weight centroids of transport zones by population and by number of buildings.

For this, we will need the census data [@INEcensus].

quarto-executable-code-5450563D

```r
#| eval: false
#| include: false
# data prep
CENSUS = st_read("/media/rosa/Dados/GIS/MQAT/original/BGRI21_LISBOA.gpkg")
CENSUS = CENSUS |> 
  select(DTMN21, SUBSECCAO, N_INDIVIDUOS, N_EDIFICIOS_CLASSICOS) |>
  rename(Mun_code = "DTMN21",
         UID = "SUBSECCAO",
         Population = "N_INDIVIDUOS",
         Buildings = "N_EDIFICIOS_CLASSICOS")
CENSUS = CENSUS |> left_join(Municipalities_names |> select(Mun_code, Municipality) |> distinct()) |> 
  select(Municipality, UID, Population, Buildings)
st_write(CENSUS, "../data/census.gpkg", delete_dsn = TRUE)
```

```{r getcensus}
#| message: false
#| fig-format: png
Census = st_read("../data/census.gpkg", quiet = TRUE)

mapview(Census |> filter(Municipality == "Lisboa"), zcol = "Population")
```

It was not that easy to estimate weighted centroids with R, as it is with GIS software.
But there is this new package [`centr`](https://ryanzomorrodi.github.io/centr) that can help us [@centr].

We need to specify the **group** we want to calculate the **mean centroids**, and the **weight variable** we want to use.

```{r centr}
library(centr)
Centroid_pop = Census |> 
  mean_center(group = "Municipality", weight = "Population")
```

We can do the same for buildings.

```{r buildings}
Centroid_build = Census |> 
  mean_center(group = "Municipality", weight = "Buildings")
```

## Compare centroids in a map

### Interactive map

```{r mapcomparecentr}
#| fig-format: png
mapview(Centroid_geo, col.region = "blue") +
  mapview(Centroid_pop, col.region = "red") +
  mapview(Centroid_build, col.region = "black") +
  mapview(Municipalities_geo, alpha.regions = 0) # polygon limits
```

See how the building, population and geometric centroids of Sintra are apart, from closer to Tagus, to the rural area.

### Static map

To produce the same map, using only `plot()` and `st_geometry()`, we need to make sure that the geometries have the same crs.

quarto-executable-code-5450563D

```r
#| eval: false
st_crs(Centroid_geo) # 4326 WGS84
st_crs(Centroid_pop) # 3763 Portugal TM06
```

So, we need to transform the geometries to the same crs.

quarto-executable-code-5450563D

```r
Centroid_pop = st_transform(Centroid_pop, crs = 4326)
Centroid_build = st_transform(Centroid_build, crs = 4326)
```

Now, to use `plot()` we incrementally add layers to the plot.

quarto-executable-code-5450563D

```r
#| fig.dpi: 300
#| fig-cap: "Static map of different centroids of Municipalities"
#| fig-format: png
# Plot the Municipalities_geo polygons first (with no fill)
plot(st_geometry(Municipalities_geo), col = NA, border = "black")

# Add the Centroids_geo points in blue
plot(st_geometry(Centroid_geo), col = "blue", pch = 16, add = TRUE) # add!

# Add the Centroid_pop points in red
plot(st_geometry(Centroid_pop), col = "red", pch = 16, add = TRUE)

# Add the Centroid_build points in black
plot(st_geometry(Centroid_build), col = "black", pch = 16, add = TRUE)
```

In the [next section](desire-lines.qmd) we will use these centroids to calculate the desire lines between them.

quarto-executable-code-5450563D

```r
#| include: false
#| eval: false
st_write(Centroid_pop, "../data/Centroid_pop.gpkg", delete_dsn = TRUE)
```
INE. 2022. “Censos 2021- XVI Recenseamento Geral da População. VI Recenseamento Geral da Habitação.” Lisboa: Instituto National de Estatística. https://censos.ine.pt/xurl/pub/65586079.
Zomorrodi, Ryan. 2024. Centr: Weighted and Unweighted Spatial Centers. https://ryanzomorrodi.github.io/centr/.